/*! \file */
/*
* Copyright (c) 2008-2011, Thomas Handorf, Nils Christian
* 
* This file is part of libScopes.
* libScopes is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* libScopes is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU Lesser General Public License
* along with libScopes.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "scopes.h"

void removeChannels(VS* where, VS* which);
void addChannels(VS* where, VS* which);
bool addChannelsNoDup(VS* where, VS* which);
bool inChannels(VS* where, VS* which);

   

/*!
\brief calculate one interation
\anchor dScopeStep
\param net the network object
\param present the currently present compounds
\param depl the currently depleted compounds
\param active array for active reactions; content is ignored, but activity of reactions is returned
\param blocking array for blocking reactions; content is ignored, but blocking of reactions is returned
\param seed compounds always present
\param block compounds always blocked (depleted)
\param chn State of Channels when network has families
\param noi bool flag to ignore inhibition

calculates one interation in the scope algorithm using the given network. 
scope, depl, seed, block are of type VC containing a value !=0 at each index whose compound is in the set and 0 otherwise. 
*/

void dScopeStep(Net* net, VC* present, VC* depl, VC* active, VC* blocking, VC* seed, VC* block, Channels* chn, bool noi){
   size_t nrr=net->numR(); //number of reactions
   size_t nrc=net->numC(); //number of compounds
   size_t i,n,n2;
   active->assign(nrr,0); // active reactions 0
   blocking->assign(nrr,0); // blocking reactions 0
   if (!noi){
      for (i=0;i<nrr;i++){ // loop for inhibited substrates
         if (!(*(net->reaMask))[i]) continue; // reaction is inactivated
         VS* coms;
         VS* ihds=net->inhibited(i); // inhibited substrates
         VS* prps=net->propagators(i); // depletion propagators
         if (ihds->size()==0) continue; // no inhibited substrates -> nothing to do
         bool fwd=true;
         bool cat=false;
         bool inh=false;
         bool prp=false;
         bool ihs;
         coms=net->substrates(i);
         for (n=0;n<coms->size();n++){
            ihs=false;
            for (n2=0;n2<ihds->size();n2++){// is current substrate an inhibited substrate?
               if ((*coms)[n]==(*ihds)[n2]) ihs=true; 
            }
            if ((!(*present)[(*coms)[n]]) || (((*depl)[(*coms)[n]]) && (!ihs))) fwd=false; // depled substrates are ok if they are the inhibited substrates
         }
         // substrate inhibition only works in forward direction -> don't need to check backward dir
         coms=net->catalysts(i);
         for (n=0;n<coms->size();n++){ // at least one catalyst must be present
            if (((*present)[(*coms)[n]]) && (!((*depl)[(*coms)[n]]))) cat=true;
         }
         if (coms->size()==0) cat=true; // if there are no catalysts or activators, cat is true
         coms=net->inhibitors(i); // inhibitors; one must be present to inhibit
         for (n=0;n<coms->size();n++){
            if (((*present)[(*coms)[n]]) && (!((*depl)[(*coms)[n]]))) inh=true;
         }
         for (n=0;n<prps->size();n++){ // if there are propagators, check whether thy are depled in order to make this reaction also blocking
            if ((*depl)[(*prps)[n]]) prp=true;
         }
         if (prps->size()==0) prp=true; // without propagators reaction is directly blocking
         if (fwd && cat && (! inh)){ // reaction is active and can inhibit
            if (prp) (*blocking)[i]=1; // blocking is dependent on propagators, activation not!!
            (*active)[i]=1; // add as active reaction here as it will not be activated in second round due to blocking -> products of the blocking reaction will be activated
         }
      }
      if (block){
         depl->assign(*block);
      } else {
         depl->assign(nrc,0);
      }
      for (i=0;i<nrr;i++){ // set newly depleted compounds
         if (!(*blocking)[i]) continue;
         VS* ihds=net->inhibited(i);
         for (n2=0;n2<ihds->size();n2++){// mark all inhibited substrates as depled
            (*depl)[(*ihds)[n2]]=1;
         }
      }
      if (chn){ // presence of chn triggers family treatment
         for (i=0;i<nrc;i++){ // propagate depletion for families
            VS* mem=net->members(i);
            if (!mem->size()) continue; // not a family
            if ((*depl)[i]){ // deplete all family members, because there exists a reaction which can deplete the whole family
               for (n=0;n<mem->size();n++){
                  (*depl)[(*mem)[n]]=1;
               }
            } else {
               for (n=0;n<mem->size();n++){
                  if ((*depl)[(*mem)[n]]) { // a family member is depleted
                     // FIXME depletion propagation will not work if only one active channel is blocked
                     // This would require a complete rework of the method based on Channels instead of compounds
                     if (!(*chn)[i]) { // family has no active channels yet
                        if ((*present)[i]){
                           (*chn)[i]=new VS(*(net->channels(i))); // copy defined channels into active channels
                        } else {
                           continue; 
                        }
                     }
                     removeChannels((*chn)[i],net->channels((*mem)[n]));
                     if ((*chn)[i]->size()==0) (*depl)[i]=1; // block family if there are no undepleted active Channels
                  }
               }
            }
         }
      }
   }
   for (i=0;i<nrr;i++){ // loop for active reactions
      if (!(*(net->reaMask))[i]) continue; // reaction is inactivated
      VS* coms;
      bool fwd=true; // substrates
      bool bck=true; // products
      bool cat=false;
      bool inh=false;
      if (net->rev(i)>=0) {
         coms=net->substrates(i);
         for (n=0;n<coms->size();n++){
            if ((!(*present)[(*coms)[n]]) || ((*depl)[(*coms)[n]])) fwd=false;
         }
      } else {
         fwd=false;
      }
      if (net->rev(i)<=0){
         coms=net->products(i);
         for (n=0;n<coms->size();n++){
            if ((!(*present)[(*coms)[n]]) || ((*depl)[(*coms)[n]])) bck=false;
         }
      } else {
         bck=false;
      }
      coms=net->catalysts(i);
      for (n=0;n<coms->size();n++){ // at least one catalyst must be present
         if (((*present)[(*coms)[n]]) && (!((*depl)[(*coms)[n]]))) cat=true;
      }
      if (coms->size()==0) cat=true; // if there are no catalysts or activators, cat is true
      if (!noi){
         coms=net->inhibitors(i); // inhibitors; one must be present to inhibit
         for (n=0;n<coms->size();n++){
            if (((*present)[(*coms)[n]]) && (!((*depl)[(*coms)[n]]))) inh=true;
         }
      }
      if ((fwd || bck) && cat && (! inh)){ // reaction is active 
         char dir=0;
         if (fwd) dir+=1;
         if (bck) dir+=2;
         (*active)[i]=dir; // add as active reaction (direction is encoded here)
      }
   }
   if (seed){
      present->assign(*seed);
   } else {
      present->assign(nrc,0);
   }
   for (i=0;i<nrr;i++){ // set newly active compounds
      if (!(*active)[i]) continue;
      VS* prod=net->products(i);
      VS* subs=net->substrates(i);
      VS* lchn=new VS();
      VS* rchn=new VS();
      VS* hchn;
      bool rfull=true; //all families on right side (products) are fully activated ((*chn)[(*prod)[n]]==NULL)
      bool lfull=true;
      if (chn){ //families?
         for (n=0;n<prod->size();n++){ // unify all product side Channels
            if ((*chn)[(*prod)[n]]){
               hchn=(*chn)[(*prod)[n]];
               if (net->members((*prod)[n])->size()){ //is family?
                  rfull=false;
               }
            } else {
               hchn=net->channels((*prod)[n]);
            }
            addChannels(rchn,hchn);
         }
         for (n=0;n<subs->size();n++){// unify all substrate side Channels
            if ((*chn)[(*subs)[n]]){
               hchn=(*chn)[(*subs)[n]];
               if (net->members((*subs)[n])->size()){ //is family?
                  lfull=false;
               }
            } else {
               hchn=net->channels((*subs)[n]);
            }
            addChannels(lchn,hchn);
         }
         for (n=0;n<prod->size();n++){ // remove non-family Channels on product side in left side channel list
                                       // note: this considers stoichiometry. If there are 2 P00001 on product side, axactly 2 P00001 are taken from substrate side
            if (!net->members((*prod)[n])->size()){ // is not a family?
               size_t j;
               for (j=0;j<net->at(i,(*prod)[n]);j++){
                  removeChannels(lchn,net->channels((*prod)[n]));
               }
            }
         }
         for (n=0;n<subs->size();n++){ // remove non-family Channels on substrate side in right side channel list
            if (!net->members((*subs)[n])->size()){ // is not a family?
               size_t j;
               for (j=0;j<(-net->at(i,(*subs)[n]));j++){
                  removeChannels(rchn,net->channels((*subs)[n]));
               }
            }
         }
         if ((*active)[i]==3){ // both directions; add family member Channels on both sides as different channels may come from both sides
            hchn=new VS();
            addChannelsNoDup(hchn,lchn);
            addChannelsNoDup(hchn,rchn);
            delete lchn;lchn=NULL;
            delete rchn;rchn=NULL;
            lchn=new VS(*hchn); // just one copy
            rchn=hchn;
            if (rfull || lfull){ // if one side is fully activated; all sides should become fully activated
               rfull=true;
               lfull=true;
            }
         }
      }
      if ((*active)[i] & 1){ //active in forward direction
         for (n=0;n<prod->size();n++){// mark all products active
            (*present)[(*prod)[n]]=1;
            if ((chn) && (net->members((*prod)[n])->size())){ // add all family member product Channels to all family products
               if (lfull){ // fully activate products
                  delete (*chn)[(*prod)[n]]; (*chn)[(*prod)[n]]=NULL;
               } else {
                  if ((*chn)[(*prod)[n]]) {
                     addChannelsNoDup((*chn)[(*prod)[n]],lchn); // add new active Channels to products (only to families)
                  } else {
                     (*chn)[(*prod)[n]]=new VS(*lchn); // copy 
                  }
               }
            }
         }
      }
      if ((*active)[i] & 2){ //active in backward direction
         for (n=0;n<subs->size();n++){// mark all substrates as active
            (*present)[(*subs)[n]]=1;
            if ((chn) && (net->members((*subs)[n])->size())){
               if (rfull){ // fully activate substrates
                  delete (*chn)[(*subs)[n]]; (*chn)[(*subs)[n]]=NULL;
               } else {
                  if ((*chn)[(*subs)[n]]) {
                     addChannelsNoDup((*chn)[(*subs)[n]],rchn); // add new active Channels to substrates (only to families)
                  } else {
                     (*chn)[(*subs)[n]]=new VS(*rchn); // copy 
                  }
               }
            }
         }
      }
      delete rchn;rchn=NULL;
      delete lchn;lchn=NULL;
   }
   if (chn){
      bool chngd=true;
      while (chngd) { // repeat while presences are changed to get propagation through multiple family hierarchy levels
        chngd=false;
        for (i=0;i<nrc;i++){ // propagate activation for families (using Channels set above)
          VS* mem=net->members(i);
          if (!mem->size()) continue; // not a family
          for (n=0;n<mem->size();n++){ // set family presence based on family member presence; this needs to be done even for already active families since new Channels may be added
              if ((*present)[(*mem)[n]]){ // family member is active
                if (!(*chn)[i]) {
                    if ((*present)[i]){
                      (*chn)[i]=net->channels(i)->copy();
                    } else {
                      (*chn)[i]=new VS();
                    }
                }
                chngd=addChannelsNoDup((*chn)[i],net->channels((*mem)[n])) || (!(*present)[i]) || chngd;
                (*present)[i]=1;
              }
          }
          if ((*present)[i]){ // activate family members
              bool nochn=(!(1 && (*chn)[i])); // no active Channels set?
              for (n=0;n<mem->size();n++){
                if (nochn || inChannels((*chn)[i],net->channels((*mem)[n]))){
                    if (!(*present)[(*mem)[n]]) chngd=true;
                    (*present)[(*mem)[n]]=1; // FIXME in principle a subset of channels would be ok and could be passed further (but we must check whether there exists a sub sub family member which is fully represented in the remaining channel set?)
                }
              }
          } else {
              if ((*chn)[i]) {
                delete (*chn)[i];
                (*chn)[i]=NULL;
              }
          }
        }
      }
   }
}

void removeChannels(VS* where, VS* which){
   size_t i,n;
   for (i=0;i<which->size();i++){
      for (n=0;n<where->size();n++){
         if ((*which)[i]==(*where)[n]){
            where->erase(n,1);
            break; // note: this removes only one element from where per one element in which
         }
      }
   }
}
bool inChannels(VS* where, VS* which){
   size_t i,n;
   size_t found=0;
   for (i=0;i<which->size();i++){
      for (n=0;n<where->size();n++){
         if ((*which)[i]==(*where)[n]){
            found+=1;
            break;
         }
      }
   }
   return found==which->size();
}
void addChannels(VS* where, VS* which){
   size_t i;
   for (i=0;i<which->size();i++){
      where->push_back((*which)[i]);
   }
}
bool addChannelsNoDup(VS* where, VS* which){ //add Channels only if a channel is not already in where
   size_t i,n;
   bool found;
   bool chngd=false;
   for (i=0;i<which->size();i++){
      found=false;
      for (n=0;n<where->size();n++){
         if ((*which)[i]==(*where)[n]){
            found=true;
            break;
         }
      }
      if (!found) {
         where->push_back((*which)[i]);
         chngd=true;
      }
   }
   return chngd;
}
/*!
\brief calculate the scope
\anchor dScope
\param net the network object
\param max maximum number of iterations
\param seed compounds always present
\param block compounds always blocked (depleted)
\param chn State of Channels when network has families

calculates the scope in the given network. 
seed, block are of type VC containing a value !=0 at each index whose compound is in the set and 0 otherwise. 
returns an array of type StateList* with the State of the network for each time step. the State contains to arrays,
State->scope (compounds present or not) State->depl (compound blocked or not).
The calculation stops if either a static or cyclic solution is detected.
*/

StateList* dScope(Net* net,size_t max,VC* seed,VC* block, Channels* chn, bool noi){
//   if ((!chn) && net->hasFamilies) chn=initChannels(net);
   if (!net->hasFamilies()) chn=NULL; // ignore channel object if there are no families
   StateList* scopes= new StateList();
   dScopeRun(net,0,max,scopes,seed,block,chn,noi);
   return scopes;
}
/*!
\brief continue calculating the scope
\anchor dScopeRun
\param net the network object
\param min minimum number of iterations
\param max maximum number of iterations
\param seed compounds always present
\param block compounds always blocked (depleted)
\param chn State of Channels when network has families

calculates the scope in the given network, continuing a previous calculation. 
May be used to change conditions at a certain time point.
scopes is the array of scopes calculated before (may be an empty array). 
seed, block are of type VC containing a value !=0 at each index whose compound is in the set and 0 otherwise. 
min forces min steps even if static or cyclic solution is encountered earlier.
set min=n and max=n to calculate exactly n steps.
*/
void dScopeRun(Net* net,size_t min,size_t max,StateList* scopes, VC* seed,VC* block, Channels* chn, bool noi){
   size_t nrc=net->numC(); //number of compounds
   VC* present;
   VC* depl;
   VC* active=new VC(nrc); // we don't track these in scopes
   VC* blocking=new VC(nrc);
   State* st;
   if (!net->hasFamilies()) chn=NULL; // ignore channel object if there are no families
   size_t size=scopes->size(); //>0 if there are already States in scopes
   scopes->cyclic=false;
   scopes->cyclelen=0;
   if (size>0){ // retrieve last State if any
      #ifdef DEBUG   
      printf("dScopeRun: found previous %d\n",(int) size);fflush(stdout); 
      #endif
      State* last=(*scopes)[size-1];
      present=last->present->copy();
      if (seed) present->OR(*seed); // add current seed
      depl=last->depleted->copy();
      if (block) depl->OR(*block);
   } else {
      present=(seed ? seed->copy() : new VC(nrc));
      depl=(block ? block->copy() : new VC(nrc));
   }
   bool finish=false;
   size_t cc=0;
   int i;
   while (((cc<max) && (!finish)) || (cc<min) ){ // step forward as long as no cyclic or static solution is reached and nr steps is less than max
      cc++;
      st=new State;
      st->present=present->copy();
      st->depleted=depl->copy();
      st->active=active->copy();
      st->blocking=blocking->copy();
      scopes->push_back(st);
      #ifdef DEBUG   
      present->print();
      depl->print();
      active->print();
      blocking->print();
      #endif
      dScopeStep(net,present,depl,active,blocking,seed,block,chn,noi);
//      present->print();
//      depl->print();
      if (!finish){ // once we now that we have a static or cyclic solution we just have to wait for cc>=min
         for (i=scopes->size()-1;i>=(int) size;i--){ // search for current State in history to detect static or cyclic solutions; only search States of current run
   //         (*scopes)[i]->scope->print();
   //         (*scopes)[i]->depleted->print();
            if ((present->isSetEq(*((*scopes)[i]->present))) && (depl->isSetEq(*((*scopes)[i]->depleted)))){
               finish=true;
               if (i<(int) scopes->size()-1){
                  #ifdef DEBUG
                  printf("dScopeRun: cyclic (%d)\n",(int) (scopes->size()-i));fflush(stdout);
                  #endif
                  scopes->cyclic=true;
                  scopes->cyclelen=scopes->size()-i;
               }
            }
         }
      }
   }
   // note: the following only works if there is no cyclic solution; that's why is commented out;
/*   while (cc<min){ // fill up the rest of min steps
      cc++;
      st=new State;
      st->scope=scope->copy();
      st->depleted=depl->copy();
      scopes->push_back(st);
   }*/
   delete active;
   delete blocking;
   delete present;
   delete depl;
}

// this gives the static-scope only
// FIXED!! we should rather get a new faster static scope function which considers catalysts and families
VC* scope(Net* net, VC* seed){
   Scope* s=sScope(net,seed);
   VC* ret=s->present;
   s->present=NULL;
   delete s; // this relies on the destructor of __scope which checks for NULL
   return ret;
}

// VC* scope(Net* net, VC* seed){
//    Channels* chn=NULL;
//    if (net->hasFamilies()) chn=new Channels(net); 
//    StateList* scopes= new StateList();
//    dScopeRun(net,0,10000,scopes,seed,NULL,chn,true); // WARNING this will at most do 10000 iterations
//    VC* ret=scopes->back()->present->copy(); //no inhibition, hence we don't need combine
//    delete chn;chn=NULL;
//    delete scopes;scopes=NULL;
//    return ret;
// }

/*!
\brief internal routine which calculates the scope
\param net the network object
\param seed an array containing the seed

internal usage. used both for scope and scopeExp
*/
template <typename T> __scope<T>* sScope(Net* net, array<T>* seed, bool exp) {
   size_t nrr=net->numR(); //number of reactions
   size_t nrc=net->numC(); //number of compounds
   size_t i,n;
   T generation=(T) 1; // the 'time' (iteration round) of the expansion process
   VS *fwdrea,*bckrea,*newcom;
   //the scope so far
   array<T>* scope=new array<T>(nrc,(T) 0);
   //the reaction sscope
   array<T>* rscope=new array<T>(nrr,(T) 0);
   //potential reactions to be added in forward direction , may contain duplicates
   fwdrea=new VS(); 
   //potential reactions to be added in backward direction , may contain duplicates
   bckrea=new VS(); 
   //new compounds from last loop
   newcom=new VS(); 
   //reserve maximum space in arrays to avoid realloc
   fwdrea->reserve(nrr);
   bckrea->reserve(nrr);
   newcom->reserve(nrc);
   newcom->clear();
   //initialize scope with seed. note: the seed array is untouched
   for (i=0;i<nrc;i++){
      if ((*seed)[i]) {
      newcom->push_back(i);
      (*scope)[i]=(T) generation;
      }
   }

   do {
      if (exp) generation++; // increase round number only if expansion is requested (type VS); otherwise it stays always 1
      fwdrea->clear();
      bckrea->clear();
      for (i=0;i<newcom->size();i++){ // for all compounds added in last round
         size_t nc=(*newcom)[i];
         VS* succs=net->succs(nc); // get succeding reactions
         for (n=0;n<succs->size();n++){ // for all succeding reactions
            size_t r=(*succs)[n];
            if ((!(*rscope)[r]) && (net->enabled(r)) && (net->rev(r)>=0)){ // is not already included, reaction is enabled in Net + check correct direction 
               fwdrea->push_back(r); //add to potential reactions
            }
         }
         //do the same for all preceeding reactions
         VS* preds=net->preds(nc); 
         for (n=0;n<preds->size();n++){
            size_t r=(*preds)[n];
            if ((!(*rscope)[r]) && (net->enabled(r)) && (net->rev(r)<=0)){ // is reaction enabled in Net + check correct direction 
               bckrea->push_back(r);
            }
         }
      }
      newcom->clear();
      for (i=0;i<fwdrea->size();i++){ //for all potential reactions that may be added in forward direction
         size_t nr=(*fwdrea)[i];
         if ((*rscope)[nr]) continue;  // already added -> next
         VS* coms=net->substrates(nr);  // get all substrates
         bool all=true; // all substrates
         for (n=0;n<coms->size();n++){  // test presence of all substrates
            if (!(*scope)[(*coms)[n]]) all=false;
         }
         bool one=false; //one catalyst
         VS* cat=net->catalysts(nr);
         if (!cat->size()) one=true; // if there is no catalyst condition is also fulfilled
         for (n=0;n<cat->size();n++){ // at least one catalyst must be present
            if ((*scope)[(*cat)[n]]) one=true;
         }
         if (all && one){
            coms=net->products(nr);
            // FIXME if generation is size_t, but T is int8, an overflow might result in 0 !??!?
            (*rscope)[nr]=(T) generation; // mark reaction as added
            for (n=0;n<coms->size();n++){ 
               size_t nc=(*coms)[n];
               if (!(*scope)[nc]) { //add products to scope and newcom if not already in scope
                  newcom->push_back(nc);
               }
            }
         }
      }
      //do the same for potential reaction that may be added in backward direction
      for (i=0;i<bckrea->size();i++){
         size_t nr=(*bckrea)[i];
         if ((*rscope)[nr]) continue;
         VS* coms=net->products(nr);
         bool all=true;
         for (n=0;n<coms->size();n++){
            if (!(*scope)[(*coms)[n]]) all=false;
         }
         bool one=false; //one catalyst
         VS* cat=net->catalysts(nr);
         if (!cat->size()) one=true; // if there is no catalyst condition is also fulfilled
         for (n=0;n<cat->size();n++){ // at least one catalyst must be present
            if ((*scope)[(*cat)[n]]) one=true;
         }
         if (all && one){
            coms=net->substrates(nr);
            (*rscope)[nr]=(T) generation; //mark reaction as added
            for (n=0;n<coms->size();n++){
               size_t nc=(*coms)[n];
               if (!(*scope)[nc]) {
                  newcom->push_back(nc);
               }
            }
         }
      }
      // FIXME we could probably do this at beginning of while loop
      for (i=0;i<newcom->size();i++){ // for all compounds added in last round
         size_t nc=(*newcom)[i];
         (*scope)[nc]=(T) generation;
      }
   } while (newcom->size()>0); // continue with next round if new compounds were added

  //clean up
  delete newcom;
  delete fwdrea;
  delete bckrea;
  
  return new __scope<T>(scope,rscope);

}
